home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / b / b.lha / B / src / bint / b1fun.c < prev    next >
C/C++ Source or Header  |  1988-11-24  |  15KB  |  624 lines
  1. /* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1985. */
  2.  
  3. /*
  4.   $Header: b1fun.c,v 1.4 85/08/22 16:48:24 timo Exp $
  5. */
  6.  
  7. /* Functions defined on numeric values. */
  8.  
  9. #include <errno.h> /* For EDOM and ERANGE */
  10.  
  11. #include "b.h"
  12. #include "b0con.h"
  13. #include "b1obj.h"
  14. #include "b1num.h"
  15.  
  16. /*
  17.  * The visible routines here implement predefined B arithmetic operators,
  18.  * taking one or two numeric values as operands, and returning a numeric
  19.  * value.
  20.  * No type checking of operands is done: this must be done by the caller.
  21.  */
  22.  
  23. typedef value (*valfun)();
  24. typedef rational (*ratfun)();
  25. typedef real (*appfun)();
  26. typedef double (*mathfun)();
  27.  
  28. /*
  29.  * For the arithmetic functions (+, -, *, /) the same action is needed:
  30.  * 1) if both operands are Integral, use function from int_* submodule;
  31.  * 2) if both are Exact, use function from rat_* submodule (after possibly
  32.  *    converting one of them from Integral to Rational);
  33.  * 3) otherwise, make both approximate and use function from app_*
  34.  *    submodule.
  35.  * The functions performing the appropriate action for each of the submodules
  36.  * are passed as parameters.
  37.  * Division is a slight exception, since i/j can be a rational.
  38.  * See `quot' below.
  39.  */
  40.  
  41. Hidden value dyop(u, v, int_fun, rat_fun, app_fun)
  42.     value u, v;
  43.     valfun int_fun;
  44.     ratfun rat_fun;
  45.     appfun app_fun;
  46. {
  47.     if (Integral(u) && Integral(v))    /* Use integral operation */
  48.         return (*int_fun)(u, v);
  49.  
  50.     if (Exact(u) && Exact(v)) {
  51.         rational u1, v1, a;
  52.  
  53.         /* Use rational operation */
  54.  
  55.         u1 = Integral(u) ? mk_rat((integer)u, int_1, 0) :
  56.                 (rational) Copy(u);
  57.         v1 = Integral(v) ? mk_rat((integer)v, int_1, 0) :
  58.                 (rational) Copy(v);
  59.         a = (*rat_fun)(u1, v1);
  60.         release((value) u1);
  61.         release((value) v1);
  62.  
  63.         if (Denominator(a) == int_1 && Roundsize(a) == 0) {
  64.             integer b = (integer) Copy(Numerator(a));
  65.             release((value) a);
  66.             return (value) b;
  67.         }
  68.  
  69.         return (value) a;
  70.     }
  71.  
  72.     /* Use approximate operation */
  73.  
  74.     {
  75.         real u1, v1, a;
  76.         u1 = Approximate(u) ? (real) Copy(u) : (real) approximate(u);
  77.         v1 = Approximate(v) ? (real) Copy(v) : (real) approximate(v);
  78.         a = (*app_fun)(u1, v1);
  79.         release((value) u1);
  80.         release((value) v1);
  81.  
  82.         return (value) a;
  83.     }
  84. }
  85.  
  86.  
  87. Hidden integer isum(u, v) integer u, v; {
  88.     return int_sum(u, v);
  89. }
  90.  
  91. Visible value sum(u, v) value u, v; {
  92.     if (IsSmallInt(u) && IsSmallInt(v))
  93.         return mk_integer(SmallIntVal(u) + SmallIntVal(v));
  94.     return dyop(u, v, (value (*)())isum, rat_sum, app_sum);
  95. }
  96.  
  97. Hidden integer idiff(u, v) integer u, v; {
  98.     return int_diff(u, v);
  99. }
  100.  
  101. Visible value diff(u, v) value u, v; {
  102.     if (IsSmallInt(u) && IsSmallInt(v))
  103.         return mk_integer(SmallIntVal(u) - SmallIntVal(v));
  104.     return dyop(u, v, (value (*)())idiff, rat_diff, app_diff);
  105. }
  106.  
  107. Visible value prod(u, v) value u, v; {
  108.     if (IsSmallInt(u) && IsSmallInt(v))
  109.         return (value)
  110.             mk_int((double)SmallIntVal(u) * (double)SmallIntVal(v));
  111.     return dyop(u, v, (value (*)())int_prod, rat_prod, app_prod);
  112. }
  113.  
  114.  
  115. /*
  116.  * We cannot use int_quot (which performs integer division with truncation).
  117.  * Here is the routine we need.
  118.  */
  119.  
  120. Hidden value xxx_quot(u, v) integer u, v; {
  121.  
  122.     if (v == int_0) {
  123.         error(MESS(400, "in i/j, j is zero"));
  124.         return (value) Copy(u);
  125.     }
  126.  
  127.     return mk_exact(u, v, 0);
  128. }
  129.  
  130. Visible value quot(u, v) value u, v; {
  131.     return dyop(u, v, xxx_quot, rat_quot, app_quot);
  132. }
  133.  
  134.  
  135. /*
  136.  * Unary minus and abs follow the same principle but with only one operand.
  137.  */
  138.  
  139. Visible value negated(u) value u; {
  140.     if (IsSmallInt(u)) return mk_integer(-SmallIntVal(u));
  141.     if (Integral(u))
  142.         return (value) int_neg((integer)u);
  143.     if (Rational(u))
  144.         return (value) rat_neg((rational)u);
  145.     return (value) app_neg((real)u);
  146. }
  147.  
  148.  
  149. Visible value absval(u) value u; {
  150.     if (Integral(u)) {
  151.         if (Msd((integer)u) < 0)
  152.             return (value) int_neg((integer)u);
  153.     } else if (Rational(u)) {
  154.         if (Msd(Numerator((rational)u)) < 0)
  155.             return (value) rat_neg((rational)u);
  156.     } else if (Approximate(u) && Frac((real)u) < 0)
  157.         return (value) app_neg((real)u);
  158.  
  159.     return Copy(u);
  160. }
  161.  
  162.  
  163. /*
  164.  * The remaining operators follow less similar paths and some of
  165.  * them contain quite subtle code.
  166.  */
  167.  
  168. Visible value mod(u, v) value u, v; {
  169.     value p, q, m, n;
  170.  
  171.     if (v == (value)int_0 ||
  172.         Rational(v) && Numerator((rational)v) == int_0 ||
  173.         Approximate(v) && Frac((real)v) == 0) {
  174.         error(MESS(401, "in x mod y, y is zero"));
  175.         return Copy(u);
  176.     }
  177.  
  178.     if (Integral(u) && Integral(v))
  179.         return (value) int_mod((integer)u, (integer)v);
  180.  
  181.     /* Compute `u-v*floor(u/v)', as in the formal definition of `u mod v'. */
  182.  
  183.     q = quot(u, v);
  184.     n = floorf(q);
  185.     release(q);
  186.     p = prod(n, v);
  187.     release(n);
  188.     m = diff(u, p);
  189.     release(p);
  190.  
  191.     return m;
  192. }
  193.  
  194.  
  195. /*
  196.  * u**v has the most special cases of all the predefined arithmetic functions.
  197.  */
  198.  
  199. Visible value power(u, v) value u, v; {
  200.     real ru, rv, rw;
  201.     if (Exact(u) && (Integral(v) ||
  202.             /* Next check catches for integers disguised as rationals: */
  203.             Rational(v) && Denominator((rational)v) == int_1)) {
  204.         rational a;
  205.         integer b = Integral(v) ? (integer)v : Numerator((rational)v);
  206.             /* Now b is really an integer. */
  207.  
  208.         u = Integral(u) ? (value) mk_rat((integer)u, int_1, 0) :
  209.                 Copy(u);
  210.             
  211.         a = rat_power((rational)u, b);
  212.         release(u);
  213.         if (Denominator(a) == int_1) { /* Make integral result */
  214.             b = (integer) Copy(Numerator(a));
  215.             release((value) a);
  216.             return (value)b;
  217.         }
  218.         return (value)a;
  219.     }
  220.  
  221.     if (Exact(v)) {
  222.         integer vn, vd;
  223.         int s;
  224.         ru = (real) approximate(u);
  225.         s = (Frac(ru) > 0) - (Frac(ru) < 0);
  226.  
  227.         if (s < 0) rv = app_neg(ru), release((value)ru), ru = rv;
  228.         if (Integral(v)) {
  229.             vn = (integer)v;
  230.             vd = int_1;
  231.         } else {
  232.             vd = Denominator((rational)v);
  233.             if (s < 0 && Even(Lsd(vd)))
  234.                 error(MESS(402, "in x**(p/q), x is negative and q is even"));
  235.             vn = Numerator((rational)v);
  236.         }
  237.         if (vn == int_0) {
  238.             release((value)ru);
  239.             return one;
  240.         }
  241.         if (s == 0 && Msd(vn) < 0) {
  242.             error(MESS(403, "in 0**y, y is negative"));
  243.             return (value) ru;
  244.         }
  245.         if (s < 0 && Even(Lsd(vn)))
  246.             s = 1;
  247.         rv = (real) approximate(v);
  248.         rw = app_power(ru, rv);
  249.         release((value)ru), release((value)rv);
  250.         if (s < 0) ru = app_neg(rw), release((value)rw), rw = ru;
  251.         return (value) rw;
  252.     }
  253.  
  254.     /* Everything else: we now know u or v is approximate */
  255.  
  256.     ru = (real) approximate(u);
  257.     if (Frac(ru) < 0) {
  258.         error(MESS(404, "in x**y, x is negative and y is not exact"));
  259.         return (value) ru;
  260.     }
  261.     rv = (real) approximate(v);
  262.     if (Frac(ru) == 0 && Frac(rv) < 0) {
  263.         error(MESS(405, "in 0**y, y is negative"));
  264.         release((value)rv);
  265.         return (value) ru;
  266.     }
  267.     rw = app_power(ru, rv);
  268.     release((value)ru), release((value)rv);
  269.     return (value) rw;
  270. }
  271.  
  272.  
  273. /*
  274.  * floor: for approximate numbers app_floor() is used;
  275.  * for integers it is a no-op; other exact numbers effectively calculate
  276.  * u - (u mod 1).
  277.  */
  278.  
  279. Visible value floorf(u) value u; {
  280.     integer quo, rem, v;
  281.     digit div;
  282.  
  283.     if (Integral(u)) return Copy(u);
  284.     if (Approximate(u)) return app_floor((real)u);
  285.  
  286.     /* It is a rational number */
  287.  
  288.     div = int_ldiv(Numerator((rational)u), Denominator((rational)u),
  289.         &quo, &rem);
  290.     if (div < 0 && rem != int_0) { /* Correction for negative noninteger */
  291.         v = int_diff(quo, int_1);
  292.         release((value) quo);
  293.         quo = v;
  294.     }
  295.     release((value) rem);
  296.     return (value) quo;
  297. }
  298.  
  299.  
  300. /*
  301.  * ceiling x is defined as -floor(-x);
  302.  * and that's how it's implemented, except for integers.
  303.  */
  304.  
  305. Visible value ceilf(u) value u; {
  306.     value v;
  307.     if (Integral(u)) return Copy(u);
  308.     u = negated(u);
  309.     v = floorf(u);
  310.     release(u);
  311.     u = negated(v);
  312.     release(v);
  313.     return u;
  314. }
  315.  
  316.  
  317. /*
  318.  * round u is defined as floor(u+0.5), which is what is done here,
  319.  * except for integers which are left unchanged.
  320.  */
  321.  
  322. Visible value round1(u) value u; {
  323.     value v, w;
  324.  
  325.     if (Integral(u)) return Copy(u);
  326.  
  327.     v = sum(u, (value)rat_half);
  328.     w = floorf(v);
  329.     release(v);
  330.  
  331.     return w;
  332. }
  333.  
  334.  
  335. /*
  336.  * u round v is defined as 10**-u * round(v*10**u).
  337.  * A complication is that u round v is always printed with exactly u digits
  338.  * after the decimal point, even if this involves trailing zeros,
  339.  * or if v is an integer.
  340.  * Consequently, the result is always kept as a rational, even if it can be
  341.  * simplified to an integer, and the size field of the rational number
  342.  * (which is made negative to distinguish it from integers, and < -1 to
  343.  * distinguish it from approximate numbers) is used to store the number of
  344.  * significant digits.
  345.  * Thus a size of -2 means a normal rational number, and a size < -2
  346.  * means a rounded number to be printed with (-2 - length) digits
  347.  * after the decimal point.  This last expression can be retrieved using
  348.  * the macro Roundsize(v) which should only be applied to Rational
  349.  * numbers.
  350.  */
  351.  
  352. Visible value round2(u, v) value u, v; {
  353.     value w, tenp;
  354.     int i;
  355.  
  356.     if (!Integral(u)) {
  357.         error(MESS(406, "in n round x, n is not an integer"));
  358.         i = 0;
  359.     } else
  360.         i = propintlet(intval(u));
  361.  
  362.     tenp = tento(i);
  363.  
  364.     v = prod(v, tenp);
  365.     w = round1(v);
  366.  
  367.     release(v);
  368.     v = quot(w, tenp);
  369.     release(w);
  370.     release(tenp);
  371.  
  372.     if (i > 0) {    /* Set number of digits to be printed */
  373.         if (propintlet(-2 - i) < -2) {
  374.             if (Rational(v))
  375.                 Length(v) = -2 - i;
  376.             else if (Integral(v)) {
  377.                 w = v;
  378.                 v = mk_exact((integer) w, int_1, i);
  379.                 release(w);
  380.             }
  381.         }
  382.     }
  383.  
  384.     return v;
  385. }
  386.  
  387.  
  388. /*
  389.  * sign u inspects the sign of either u, u's numerator or u's fractional part.
  390.  */
  391.  
  392. Visible value signum(u) value u; {
  393.     int s;
  394.  
  395.     if (Exact(u)) {
  396.         if (Rational(u))
  397.             u = (value) Numerator((rational)u);
  398.         s = u==(value)int_0 ? 0 : Msd((integer)u) < 0 ? -1 : 1;
  399.     } else
  400.         s = Frac((real)u) > 0 ? 1 : Frac((real)u) < 0 ? -1 : 0;
  401.  
  402.     return MkSmallInt(s);
  403. }
  404.  
  405.  
  406. /*
  407.  * ~u makes an approximate number of any numerical value.
  408.  */
  409.  
  410. Visible value approximate(u) value u; {
  411.     double x, e, expo;
  412.     register int i;
  413.     if (Approximate(u)) return Copy(u);
  414.     if (IsSmallInt(u))
  415.         x = SmallIntVal(u), expo = 0;
  416.     else if (Rational(u)) {
  417.         rational r = (rational) u;
  418.         integer p = Numerator(r), q;
  419.         double xp, xq;
  420.         int lp, lq;
  421.         if (p == int_0)
  422.             return Copy((value)app_0);
  423.         q = Denominator(r);
  424.         lp = IsSmallInt(p) ? 1 : Length(p);
  425.         lq = IsSmallInt(q) ? 1 : Length(q);
  426.         xp = 0, xq = 0;
  427.         expo = lp - lq;
  428.         if (IsSmallInt(p)) { xp = SmallIntVal(p); --expo; }
  429.         else {
  430.             for (i = Length(p)-1; i >= 0; --i) {
  431.                 xp = xp*BASE + Digit(p, i);
  432.                 --expo;
  433.                 if (xp < -Maxreal/BASE || xp > Maxreal/BASE)
  434.                     break;
  435.             }
  436.         }
  437.         if (IsSmallInt(q)) { xq = SmallIntVal(q); ++expo; }
  438.         else {
  439.             for (i = Length(q)-1; i >= 0; --i) {
  440.                 xq = xq*BASE + Digit(q, i);
  441.                 ++expo;
  442.                 if (xq > Maxreal/BASE)
  443.                     break;
  444.             }
  445.         }
  446.         x = xp/xq;
  447.     }
  448.     else if (Integral(u)) {
  449.         integer p = (integer) u;
  450.         if (Length(p) == 0)
  451.             return Copy((value)app_0);
  452.         x = 0;
  453.         expo = Length(p);
  454.         for (i = Length(p)-1; i >= 0; --i) {
  455.             x = x*BASE + Digit(p, i);
  456.             --expo;
  457.             if (x < -Maxreal/BASE || x > Maxreal/BASE)
  458.                 break;
  459.         }
  460.     }
  461.     while (expo < 0 && fabs(x) > BASE/Maxreal) ++expo, x /= BASE;
  462.     while (expo > 0 && fabs(x) < Maxreal/BASE) --expo, x *= BASE;
  463.     if (expo != 0) {
  464.         expo = expo*twologBASE + 1;
  465.         e = floor(expo);
  466.         x *= .5 * exp((expo-e) * logtwo);
  467.     }
  468.     else
  469.         e = 0;
  470.     return (value) mk_approx(x, e);
  471. }
  472.  
  473.  
  474. /*
  475.  * numerator v returns the numerator of v, whenever v is an exact number.
  476.  * For integers, that is v itself.
  477.  */
  478.  
  479. Visible value numerator(v) value v; {
  480.     if (!Exact(v)) {
  481.         error(MESS(407, "*/ on approximate number"));
  482.         return zero;
  483.     }
  484.  
  485.     if (Integral(v)) return Copy(v);
  486.  
  487.     return (value) Copy(Numerator((rational)v));
  488. }
  489.  
  490.  
  491. /*
  492.  * /*v returns the denominator of v, whenever v is an exact number.
  493.  * For integers, that is 1.
  494.  */
  495.  
  496. Visible value denominator(v) value v; {
  497.     if (!Exact(v)) {
  498.         error(MESS(408, "/* on approximate number"));
  499.         return zero;
  500.     }
  501.  
  502.     if (Integral(v)) return one;
  503.  
  504.     return (value) Copy(Denominator((rational)v));
  505. }
  506.  
  507.  
  508. /*
  509.  * u root v is defined as v**(1/u), where u is usually but need not be
  510.  * an integer.
  511.  */
  512.  
  513. Visible value root2(u, v) value u, v; {
  514.     if (u == (value)int_0 ||
  515.         Rational(u) && Numerator((rational)u) == int_0 ||
  516.         Approximate(u) && Frac((real)u) == 0) {
  517.         error(MESS(409, "in n root x, n is zero"));
  518.         v = Copy(v);
  519.     } else {
  520.         u = quot((value)int_1, u);
  521.         v = power(v, u);
  522.         release(u);
  523.     }
  524.  
  525.     return v;
  526. }
  527.  
  528. /* root x is computed more exactly than n root x, by doing
  529.    one iteration step extra.  This ~guarantees root(n**2) = n. */
  530.  
  531. Visible value root1(v) value v; {
  532.     value r = root2((value)int_2, v);
  533.     value v_over_r, theirsum, result;
  534.     if (Approximate(r) && Frac((real)r) == 0.0) return (value)r;
  535.     v_over_r = quot(v, r);
  536.     theirsum = sum(r, v_over_r), release(r), release(v_over_r);
  537.     result = quot(theirsum, (value)int_2), release(theirsum);
  538.     return result;
  539. }
  540.  
  541. /* The rest of the mathematical functions */
  542.  
  543. Visible value pi() { return (value) mk_approx(3.141592653589793238463, 0.0); }
  544. Visible value e() { return (value) mk_approx(2.718281828459045235360, 0.0); }
  545.  
  546. Hidden value trig(v, fun, zeroflag) value v; mathfun fun; bool zeroflag; {
  547.     real w = (real) approximate(v);
  548.     double expo = Expo(w), frac = Frac(w), x, result;
  549.     extern int errno;
  550.     if (expo <= Minexpo/2) {
  551.         if (zeroflag) return (value) w; /* sin small x = x, etc. */
  552.         frac = 0, expo = 0;
  553.     }
  554.     release((value)w);
  555.     if (expo > Maxexpo) errno = EDOM;
  556.     else {
  557.         x = ldexp(frac, (int)expo);
  558.         if (x >= Maxtrig || x <= -Maxtrig) errno = EDOM;
  559.         else {
  560.             errno = 0;
  561.             result = (*fun)(x);
  562.         }
  563.     }
  564.     if (errno != 0) {
  565.         if (errno == ERANGE)
  566.             error(MESS(410, "the result is too large to be representable"));
  567.         else if (errno == EDOM)
  568.             error(MESS(411, "x is too large to compute a meaningful result"));
  569.         else error(MESS(412, "math library error"));
  570.         return Copy((value)app_0);
  571.     }
  572.     return (value) mk_approx(result, 0.0);
  573. }
  574.  
  575. Visible value sin1(v) value v; { return trig(v, sin, Yes); }
  576. Visible value cos1(v) value v; { return trig(v, cos, No); }
  577. Visible value tan1(v) value v; { return trig(v, tan, Yes); }
  578.  
  579. Visible value atn1(v) value v; {
  580.     real w = (real) approximate(v);
  581.     double expo = Expo(w), frac = Frac(w);
  582.     if (expo <= Minexpo + 2) return (value) w; /* atan of small x = x */
  583.     release((value)w);
  584.     if (expo > Maxexpo) expo = Maxexpo;
  585.     return (value) mk_approx(atan(ldexp(frac, (int)expo)), 0.0);
  586. }
  587.  
  588. Visible value exp1(v) value v; {
  589.     real w = (real) approximate(v);
  590.     real x = app_exp(w);
  591.     release((value)w);
  592.     return (value) x;
  593. }
  594.  
  595. Visible value log1(v) value v; {
  596.     real w = (real) approximate(v);
  597.     real x = app_log(w);
  598.     release((value)w);
  599.     return (value) x;
  600. }
  601.  
  602. Visible value log2(u, v) value u, v;{
  603.     value w;
  604.     u = log1(u);
  605.     v = log1(v);
  606.     w = quot(v, u);
  607.     release(u), release(v);
  608.     return w;
  609. }
  610.  
  611. Visible value atn2(u, v) value u, v; {
  612.     real ru = (real) approximate(u), rv = (real) approximate(v);
  613.     double uexpo = Expo(ru), ufrac = Frac(ru);
  614.     double vexpo = Expo(rv), vfrac = Frac(rv);
  615.     release((value)ru), release((value)rv);
  616.     if (uexpo > Maxexpo) uexpo = Maxexpo;
  617.     if (vexpo > Maxexpo) vexpo = Maxexpo;
  618.     return (value) mk_approx(
  619.         atan2(
  620.             vexpo < Minexpo ? 0.0 : ldexp(vfrac, (int)vexpo),
  621.             uexpo < Minexpo ? 0.0 : ldexp(ufrac, (int)uexpo)),
  622.         0.0);
  623. }
  624.